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Figure 1: a. Vertical response functions of 10-, 20-, 30-, 60-, 
and 90-in. median curves at 0 background conductivity, b. 60- 
in. vertical response function for all background conductivities. 



same verticaJ resolution; median depths of investigation 
are 10, 20, 30, 60, and 90 in. The weight coefficients that 
produce these logs are a function of background conduc- 
tivity and are designed to maintain constant radial and 
vertical resolution over a wide range of background con- 
ductivity. The logs have a tight vertical-invariant focus 
into the formation independent of formation conductiv- 
ity. The processed logs have a constant radial median 
response over the same range of background conductivity. 
For example, a 30-in. log is produced by formation eddy 
currents, half of which are inside and half outside of 30 
in. The set of five processed log curves is closely vertically 
matched over narrow depth intervals around vertical step 
clianges in conductivity. The uncertainties associated with 
vertical bsd-boundary effects are minimized. This effec- 
tively deconvolves the vertical dimension from the prob- 
lem except in the immediate vicinity of bed boundaries. 



Thus, by design— and prior removal of borehole fields— 
when the log curves separate, invasion has occurred. 

The five processed logs are inverted for the underlying 
radial conductivity profile. The relationship between the 
processed logs and the radial conductivity profile is a 
discretized set of pseudo- radial geometric factors stored 
on a computer for a range of background conductivities. 
The radial profile is a four-parameter Butterworth func- 
tion, which includes a transition region of resistivity with 
arbitrary inner and outer radii. The four-parameter pro- 
file at each depth is obtained using a modified Levenberg- 
Marquardt (LM) nonlinear iterative solver. The LM al- 
gorithm minimizes the ,weighteo^su^ 

nro^emjp&u^ 

The minimum leasts defihes the four param- 

eters and hence the radial profile. 

An example set of response functions ( 1 00 vertical sam- 
ples) that are linear combinations of all AIT arrays is 
shown in Figures la and lb. Figure la shows that all 
curves 10, 20, 30, 60, and 90 in. are designed to have the 
same vertical response. Figure lb is the vertical response 
of a 60-in. median filter set for all background conduc- 
tivities. These two figures demonstrate the degree to which 
it is possible to focus a fixed vertical response independent 
of background conductivity. Although not shown, me- 
dian radial depth of investigation is also maintained. The 
ability to achieve this invariant and tight focus of the 
synthesized response is the justification for the decoupling 
of the radial processing. Details of median curve-filter 
design are available (Hunka et al., 1990). 

PSEUDO-RESPONSE FUNCTIONS 

From a measurement point of view, the apparent con- 
ductivity of the rth median curve, is related through 
response function theory to the underlying formation by 



(1) 



where a m is the digitized value of the radial conductivity 
profile at radius p = p mt N R is the number of radial pixels, 
which is taken to be between 60 and 100 with a 2-in. 
radial increment, and the weighting function w tm is to be 
determined. 

Equation (1) as it stands is not useful for the estimation 
problem. The complete definitions are lacking, the phys- 
ics of invasion that determines the shape of the radial 
profile is not imposed, and the tool physics that relates 
the response to the formation is missing. The tool physics 
enters through pseudo-response functions based upon rig- 
orous radial field expansions of Maxwell's equations (Wait, 
1984). The concept of pseudo-response functions can be 
found in Tinman (1986). For this application, note that 
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■adial conductivity profile containing a single step change 
,u conductivity, Atr, i.e., 

oip) - cr 0 + Acrw(p - Po). ( 2 ) 
defines in a natural way the pseudo-response function 
J&(p) for the yth induction array 0' - U 2> • • • » N c) through 
the apparent-conductivity relation 



ffL + A<r f JKpMp - Po) dp- 
Jo 



(3) 



Here is the apparent conductivity for formation con- 
ductivity a 0 . Noting that the derivative of the unit step 
function n(p - p 0 ) is the impulse function 5(p - p 0 ), it 
follows from equation (3) that 



1 da I 



(4) 



For accurate computation of «r„ it is advantageous to 
define a rodial cumulative pseudo-response function 

J&Po)** 

4(Po)= I -MM*- 



(5) 



The advantage of this response function definition is that 
the area from a fixed radius to infinity is computed ac- 
curately as the primary numerical result. Beyond the in- 
vaded zone, this is the required quantity to yield optimum 
<x, estimates. For numerical computation with error of 
order (.r al - <r fl2 ) 2 , a central-difference scheme is used. 
Thus it is necessary to compute two apparent conductiv- 
ities, a ol and a a2 , viz.: 

* ii = ° l + X" jKf>w% " ^ *• (6a) 

<fL - + J •WG»)[»i + AffM t" - A>) ~ ff ol (6b) 

The difference of the two equations (6a) and (6b) defines 
the desired partial cumulative response, i.e., 
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Figure 2: Partial cumulative radial response function A7(p) as 
defined by equation (8) for deep reading array and background 
conductivities, as shown in legend- 
Note that iV - -/J(0). An example normalized partial- 
cumulative radial-response function computed by equa- 
tion (8) for a deep-reading array is given in Figure 2. In 
the next section, the pseudo-response function will be 
used to deconvolve the formation conductivity. In the 
appendix, analytic forms of the normalization function 
are given for both the dipole and more general loop cases. 

RADIAL DECONVOLUTION 
In this section, it is assumed that the formation resis- 
tivity varies only radially. Using the response functions 
defined in the previous section, the apparent conductivity 
measured by the jth channel is simply 



(7) 



a j= £ K'(p)*(p) dp. (10) 

For numerical purposes, it is convenient to define the 
radial conductivity profile a(p) as the finite-element 
model 



(H) 



In equation (11), A m (p) is the max triangular finite element 
defined as 



Equation (7) is the numerical definition of the pseudo- 
response function, which is used to compute it given a 
forward model. Skin correction is accomplished by nor- 
malizing response functions to have total unity gain. The 
definition is 

KJ(p) = -P0>)/A»(<ro), (»> 

where 

NK ** = J7 j,{p) dp - (9) 



A m (p) 



[l ~ IP- 

\i - !p- 
o, 



- p. 
p< 



,|/A«, 



if p m _, < p < p«; 

if Pm < P < Pm+U 

otherwise. 

(12) 



where m - 1, 2, . . . , N R - 1, A m = Pm - Pm-u and P? - 
0. Note in equation ( 1 0) that because of the normalization 
of the kernel functions, A?, it is not necessary to include 
an explicit additive background conductivity. Note also 
lhat the /> axis subdivision need not be uniform, it is only 
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necessary that it be ordered, p 0 < p, < . . . < p„ R . The 
first and last triangular functions are special cases given by 



__ , . - [p - p 0 |/A„ if p 0 < p < p,; 
~ in otherwise. 



I* 



(13) 



- Ip - p^I/a^ if p^. ( < p < P/v„; 

otherwise. 

(14) 



The linear combination of the triangular functions in 
equation (1 1) yields a piecewise linear function over the 
entire finite-element interval (p 0 , p NR ) and hence yields is 
a higher-order expansion than one with piecewise con- 
stant elements. Beyond p HR , the unit step function in equa- 
tion (11) forces the profile to the constant p NR . Substitu- 
tion of equation (11) into equation (10) yields . 



where 



F m ~ K m {p)KKp) dp. 



(15) 
(16) 



The finite-element matrix /* m can be written in terms of 
partial radial-response functions defined by equation (7), 
i.e., let 

G> m «* ^(Pj/^o). (17) 

Then the finite-element matrix P* m has the explicit rep- 
resentation 



1 m 



(1 - ^.)/2> 
(CV_, - GU,)/2, 
(Gj^., + G&J/2, 



if m « 0; 

if 0 < m < A^; (18) 
if m « iV*. 



The partial cumulative response functions G* m (<r 0 ) are 
computed and stored in table form on computer disk. 
The table consists of uniform 2-in. steps from 2 to 200 
in. (m = 1,2,..., 100) for each array (f » 1, 2, . . . , 28) 
and background conductivities cr 0 558 0, 0.01,-0.02, 0.05, 
0.1, 0.2. 0.5, I, 2, 3.5, and 5 S/m. A check on the nor- 
malization is made by choosing <r m = a c for all m. Then 
all curves should read the correct homogeneous value c c . 

I** 

= + <%)/2, 

(19) 

In citation (19), the fact Go — 1, which follows from 
equation (8), is needed. As an example of the response 
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Figure 3: Real part of differential radial response functions for 
ail arrays defined by finite element matrix P J m (cr 0 ) for <r 0 » 1 S/ 
m (see equation [16]). The response function is based upon 
triangular finite-element shape functions A m (p). 

functions, P J mJ the matrix is shown in Figure 3 for back- 
ground of 1 S/m for all arrays j. The legend labels raw 
array responses 1-8 in order of increasing depth of in- 
vestigation. The abscissa is m converted to inches. An 
important point illustrated by Figure 3 is that it is more 
difficult to focus the response radially beyond about 30- 
50 in. than it is to focus the response vertically because 
the focus is achieved in the static limit of Maxwell's equa- 
tions rather than by interference effects that are respon- 
sible for high resolution in microwave imaging where the 
size of the narrow beam-forming antennas is more than 
1 wavelength. Microwave imaging intrinsically allows 
much higher resolution as the distance from the source 
is increased. In high-loss media, as encountered in well 
logging, skin-effect is not negligible, and it is necessary to 
use induction methods. Then the ability to focus the qua- 
si-static field is based upon the much more difficult and 
limited ability to shape near-field source singularities. 
This ability degrades rapidly as one moves away from 
the borehole, as is evident in Figure 3. 

As discussed in the introduction, the primary mea- 
surements are combined to produce focused apparent- 
conductivity median logs. Median logs in the case of thick 
beds are produced by linear combinations of the primary 
measurements of the form 



He 



(20) 



The index j is reserved for actual sonde data, and the 
index / is associated with the processed or median curves. 
The. weights &j(<r 0 ) produce a focused log with median 
depths of investigation of 10, 20, 30, 60, and 90 in., as i 
= 1, 2, 3, 4, and 5. The weights are designed to have 
identical step responses for all curves (i) and background 
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Figare 4: Pour.para.eter maximally flar profiles for M - 2. 
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Figure 55: Four-parameter maximally fiat profiles for N = 4. 
r, - w(l - 2/AO, 

r 2 - w(l + 2/N). (25) 

Radii r t and r 2 are defined as the intersections of straight 
lines atssociated with, respectively, the inner and outer 
horizontal asymptotes, and the line passing through the 
midpoint (w, [<r t - <r w ]/2) and tangent to the profile at 
the same point. The order parameter M which is pro- 
portional to the slope of the transition, and the midpoint 
transition radius w are nonlinear parameters, whereas the 
asymptotes a m and cr, are linear parameters of the max- 
imally flat model. 

The parametric profile can easily be extended to allow 
for an annulus. As an example, a six-parameter profile 
<r 6 (p) witn toree conductivity parameters c^, <r flM , o, f two 
transition radii, w t , w 2 (where 0 < w t < Wj), and one 
slope parameter N is 

a 6 (p) = a t + (a m - 0/lP> "i> W 

+ iff an ~ V<)Af>> W 2> ( 26 > 

where ./(p, w. N) = 1/[1 + (p/"/*]. 
3 
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The limit of large N is an important limiting case that 
allows the identification of the significance of the three 
conductivities appearing in equation (26). Then equation 
(26) reduces to a rectangular step profile, i.e., 

if 0 < p < w,; 
if w t < p < w 2 ; (27) 
if p > w 2 . 



lim a 6 (p) — 4 <r„ 
I*.. 



An example <r 6 (p) profile with w, = 20 and w 2 = 30 is 
given in Figure 7 for several values of N as indicated in 
the legend. The four-parameter method developed in this 
paper could apply equally well to profiles given by equa- 
tions (24) and (26) or any number of radial regions. The 
fundamental limitation is the ability of induction data to 
resolve more than perhaps six parameters, 

NONLINEAR PARAMETER ESTIMATION 

The parametric conductivity profile model is nonlinear 
in the transition radius w and the order N. Thus a non- 
linear estimation method is required. The method of 
choice, resulting from work with several methods, is a 
modified LM algorithm adapted from Numerical recipes 
(Press et al., 1986). The idea is to match the apparent- 
conductivity log data <r ai with maximally flat model-based 
estimates £*/(ni) (see equation [22]) for the / - 1,2,..., 
M c by adjusting the formation conductivity profile. The 
hat notation of statistical estimation theory denotes an 
estimate. The argument m of the estimate is the estimated 
model vector whose components in the case of a four- 
parameter profile are given by 

m 2 = v f > 
m y = vv, 

m 4 - N. (28) 
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Figure 6: Four-parameter maximally flat profiles for W= 16. 
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Figure 7: Six -parameter maximally flat profiles with conduc- 
tive annulus with break points of w, = 20 and w 2 « 30 in. 
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>e estimat ion method does not depend upon the specific 
.mber of measurements, M„ or model parameters, N p . 
Of course, how well the method performs does depend 
on M c and ,V P . The model parameters enter through the 
radial pixel estimates b m of the formation conductivity. 
From equation (23), the relation is 



/n— O 



(29) 



these deficiencies by adding a numerical parameter X, 
which improves the conditioning problem. It also enables 
the iterative procedure to alternate between a gradient 
search when the current model estimate m is far from the 
minimum, and therefore approximation (31) is not ap- 
plicable, and a quadratic convergence solution approach- 
ing equation (34) near the minimum. The LM algorithm 
itcratively solves the modified system 

. = nw+ L-l-Wny)], (35) 



where b m ~ o Np {f> m )- 

In equation (29), <x Sp <j>) is an N p parameter maximally 
flat profile such as equations (24) or (26), and p m is a 
member of the discretization of the p axis as defined 
between equations (12) and (13). In the usual application 
of the LM algorithm, and in terms of these definitions, 
the model vector m is the vector that minimizes the x 2 
fit to the da ta. The definition is 



X 2 (m) - 2 I fa- - *JpW«A\ 



(30) 



where a ol is. the ith input median apparent-conductivity 
data channel. The conductivity <r, is the associated stan- 
dard deviation of the ith data channel and is an a-priori 
estimate of the reliability of the data; alternatively r« can 
epresen t the ability of the underlying model for the radial 
conductivity profile to fit the actual profile. 

The method requires knowledge of the gradient (first 
derivative) and the second partial derivative matrix (Hes- 
sian matrix) of with respect to the model vector m. 
The goodn ess of fit x*(m) is approximated by a quadratic 
surface in the model parameters m k 



1 



(31) 



m„ 



m„ 



where 



(36) 



A necessary condition for the minimum rn^m, of x 2 (m) 
is that the gradient with respect to the model parameters 
be zero: 

V X 2 (m) - 0, (32) 

which from equation (31) yields the condition 

Hm,~p. (33) 

An iterative equation for the minimum is found by setting 
the gradient of equation (3 1 ) to zero and then substituting 
£ from equation (33) into the result to obtain the iterative 
solution 

(34) 



and where 5^ is the Kronecker delta, which is 1 if n = 
m and zero otherwise. When X is small, equation (35) 
reverts to equation (34), which is appropriate near the 
minimum. However when X is large, the matrix L and its 
inverse become diagonal. In this limit, equation (35) be- 
comes a gradient method that chooses a new estimate 
model vector, m new , to be downhill from the previous 
estimate, nw The algorithm contains a loop that solves 
equation (35). After each iteration, x^ncJ i* compared 
with xHm^). If it is smaller, X is decreased by a factor of 
10, otherwise it increases X by a factor of 10, Then 
is set equal to and iteration continues until con- 
vergence is achieved. A more explicit form of equation 
(34) is found by computing the gradient and Hessian ma- 
trix from equation (30). The results are 

- 20*, 

~f <r? dm k dm, 9 



~ 2H„. (37) 
Define the change in the model vector as 

5m = m new - nw (38) 
Then the matrix equation (35) has the more explicit form 

(39) 



f-1 



As it stands, equation (34) is not suitable. Near the 
minimum it has quadratic convergence but may fail to 
find the minimum unless the initial estimate is close 
enough. In addition, the Hessian matrix H can be singular 
or ill conditioned. The LM algorithm cleverly overcomes 



For cases of deep invasion (w > 30), the induction 
measurement is not sensitive to the shape of the forma- 
tion conductivity profile controlled by model parameter 
N as can be determined by examining the response func- 
tions defined by equation (5). Another poorly resolved 
case occurs in tight rocks, where the radial profile can be 
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featureless. In such cases, which often occur, matrix equa- 
tion (39) is poorly conditioned and sometimes singular. 
In least-square estimation methods, the system matrix is 
often modified to improve its condition by adding a con- 
straint matrix such as minimum energy or curvature. 
Twomey (1977) gives several examples showing the so- 
lution enhancement possible by conditioning. If the sys- 
tem is linear, the optimum least-square method with min- 
imum norm solution is obtained by singular value 
decomposition. 

For linear systems, certain classes of constraint yield a 
unique solution to the model vector m. In the case of 
minimum energy, uniqueness can be proven (Howard et 
al., 1983). Uniqueness occurs because the conditioned 
system matrix can be shown to be nonsingular. As a fur- 
ther benefit, conditioned methods are a means to apply 
additional physical constraints on the model vector m. 
The LM problem is nonlinear and has a special x 2 struc- 
ture, so dynamically constrained LM (DCLM) is an ex- 
tension ofthe method. For the induction problem, a sim- 
ple but efFective constraint modification of LM results 
from appending a dummy data channel to the x 2 variable, 



i.e., 



X 2 (m) - xBXm) + C(m), 



(40) 



where the constraint C(m) is chosen to be a simple qua- 
dratic form in the model parameters 

C(m) = |[<w e +, - ^ c+ j(m)]/o-^ +l | 2 , (41) 

where 

The constraint tends to keep the four parameters (see 
equation [28]) in reasonable ranges 

-Vmk + Mom < m k < ow + m Qkt (42) 

for k = 1,2,..., N pt An important additional hard con- 
straint for this profile estimation problem is that all pa- 
rameters are positive. Because the LM method is nonlin- 
ear, positi vity is guaranteed by the change of model 
variables 



which modifies matrix equation (39): 



2 U k Aq, = P\ 



(43) 



(44) 



where 



L'*, = AqtfJLl + MI„). 



and where H w and ft k are given by equation (37). The 
relative emphasis ofthe solution on the constraint C(m) 
is determined by <r Nc + % . It is not possible to predict (within 
an order of magnitude) the model data fit Xo(m) achieved 
on an arbitrary input data set. This observation motivates 
the dynamic adjustment of u aNc ^ % The dynamic adjust- 
ment of the constraint weight attempts to keep each ma- 
trix equation in the nonlinear iteration method well 
conditioned. In the process of convergence, the data-de- 
pendent term Xo(™) gradually, iteration by iteration, be- 
comes smaller by several orders of magnitude. For a good 
solution during the iteration, it is necessary for the con- 
ditioning to become smaller at the same rate. This is an 
additional complication compared to linear conditioned 
estimation. The DCLM algorithm also deals with cases 
where the model may not fit the data and smoothly tracks 
a minimum x 2 solution into poorly conditioned param- 
eter regimes. 

Dynamic conditioning results by choosing cr^, such 
that the relative weight of the constraint to model 
data fit xil(m) ratio is kept at a constant value y 

C(m)/x2(m) = y. (45) 

The DCLM iteration path to the minimum x 2 solution 
can go through regions of large condition number c N of 
the matrix L' w denned by equation (44). In general too, 
the problem is ill conditioned, so that the choice of matrix 
algorithm for solving equation (44) is critical. Good ma- 
trix algorithms for this purpose are Unpack (Dongarra et 
al., 1979) double-precision symmetric indefinite factor- 
ization, DSICO, and solver, DSISL (which also report the 
condition number). 

In summary, the new DCLM-algorithm iteration scheme 
is 1) initialize parameter vector m; 2) solve equation (44) 
for the new model parameters m^; 3) determine the 
relative-constraint weighting numerical parameter o> +l 
according to equation (45), where o N<f¥l is defined by equa- 
tion (41) (in the numerical results to follow, the ratio 
parameter y is chosen to be and 4) check results to 
determine if the solution has converged, and if not, repeat 
the iteration. 



NUMERICAL RESULTS 

Several cases of invasion profiles in thick beds are shown 
in Figures 8-15. In each of these, the upper half of the 
figure is the result plotted in log resistivity flog I0 {i?]), 
where R = \/a is in ohm-m. The same data are plotted 
in the lower half of the figure with a linear conductivity 
a scale in mS/m. The abscissa is the radial coordinate r 
between the borehole radius and 1 20 in. Ofthe two equiv- 
alent representations, the conductivity is actually com- 
puted and more closely reflects the sensitivity ofthe mea- 
surement. 
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Figure 8: Thick-bed simulation where large . dot* are ^ median 
curves, cloned line is model formation, dashed bne* mjl 
estimate, and solid line is DCLM estimate for x2 - I -09 x 10°, 
7.56 >: 10*, N - 16.1, w = 9.5, and iterations - 59. 
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Figure 10: Thick-bed simulation where large dote are median 
cu^et dotted line is model formation, dashed Ime is ^miud 
estimate, and solid line is DCLM estimate for xJ- 2.08 x 
c% 2.25 x N = 6.0, w - 33.2, and interactions - 

24. 
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Figure 9: Thick-bed simulation where large dots are median 
curves, dotted line is model formation, dashed unejs mUial 
estimate and solid line is DCLM estimate for x? - 4.64 x 
?0™£ - 5.85 x 1(P, AT = 3.69, - - 19.9. and iterations - 
26. 



r, (in.) r 2 (in.) <r, (mS/m) <*„ (mS/m) 



60 

r Cncheal 

Figure 11- Thick-bed simulation where large dots are median 
cuSes, dotted Une is model formation, dashed line is initial 
estimate, and solid line is DCLM estimate for x2 - L60 x 
?0™X= M3 x 10 4 , N ■= 15.6, 19.6, and interactions = 
47. 
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Figure 12: Thick-bed simulation where large dots arc median 
curves, dotted line is model formation, dashed line is initial 
estimate, and solid line is DCLM estimate for Xo "= 3.45 x 
10°, iv- 6.17 x 10 $ , AT « 1.68, w=22.0, and iterations - 41. 
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Figure 13: Thick-bed simulation where large dots are median 
curves, dotted line is model formation, dashed line is initial 
estimate, and solid line is DCLM estimate for Xo ^ 1.61 * 
I0" 1 , c N - 7.54 x I0 3 , N= 6.6, w= 18.6, and iterations = 22. 
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The five distinct open dots are the 10-, 20-, 30-, 60-, 
and 90-in. median-log input values. The dotted line is 
the input profile, the dashed line is the initial estimate, 
which is simply a four-parameter curve fit to the four 
median points, and the solid line is the DCLM estimate. 
These figures also give in table form the four parameters 
from equation (28), and both the input and estimated 
values of the two characteristic radii r, and r 2 defined by 
equation (25), a xo , and a t . The captions also include the 
Xo fit and condition number c„ of the DCLM matrix U 
of equation (44). The quality of the solution is provided 
for the model parameters at each depth. Good solutions 
have small condition numbers^ less than perhaps 10 4 , 
indicating numerical reliability. How well the model fits 
the data is more important and is measured by Xo from 
equation (40). When Xo < 2, the fit is good. 

In several instances, Figures 8-1 5 show that the median 
logs, although indicative of the shape of the profile, are 
not necessarily quantitatively correct because they are not 
point estimates but rather are center-weighted averages 
as their descriptor "median** suggests. Therefore, radial 
profiling schemes that do little more than curve fit to the 
median log values (as does the displayed dashed curve 
used for initializing DCLM) can have unacceptably large 
error. 

One application of the shape of the invasion profile 
obtained with DCLM, which uses r,, r 2 , a xo , and a ty is to 
estimate the volume of mud filtrate per unit depth that 
has displaced the connate fluids. Figure 16 is a three- 
dimensional perspective of a simulated log through a fi- 
nitely thick 20-ft bed with conductive invasion. The sim- 
ulated data are computed by the semianalytic method 
(Chew et al., 1984). The exact model and estimated typ- 
ical values for <r„ r„ and r 2 are also included in table 
form for each of the three beds. Notice that in many of 
the shallow invasion results, the accuracy of the a, esti- 
mate is better than a^, The method is designed to have 
better <r, accuracy because the maximally flat-parameter 
model improves the sensitivity of the measurement to <r,. 
The relatively less accurate <r xo for shallow invasion cases 
is related to the radial response functions, which are small 
near the borehole (as can be seen from Figure 3). None- 
theless, the parametric processing improves the induc- 
tion-based Rjb, measurement. 

No inversion study is complete without characteriza- 
tion of noise amplification, particularly when the pro- 
cessing is nonlinear and the linear steps in the iterations 
can have large condition numbers. Noise is added to four 
input curves in the following way. Noise of 10, 5, 2, and 
i mS/m standard deviation is added to the 1 0-, 20-, 30-, 
and 60-in. curves. The noise is computed by a Box-Muller 
transformation of uniformly distributed white noise. This 
transformation produces normally distributed random 
noise (Press et al., 1986, pp. 195-203). A large number, 
N„ t of such noise samples are added to the input curves, 
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Figure 14: Thick-bed simulation where large dots are median 
curves, dotted line is model formation, dashed Line is initial 
estimate, ar.d solid line is DCLM estimate for y& - 3.4 x 
10- 1 c*«2.24x 10 J , N=> 3.6. 37.3, and iterations - 151. 
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Figure 15: Thick-bed simulation where large dots are median 
curves, dotted line is model formation, dashed line is initial 
estimate, and solid line is DCLM estimate for xi - 3.47 x 
l0 -i f Cn = i.76 x 10*, N = 6.5, w = 34.8, and iterations « 24. 



r» (in.) 


r 2 (in.) 


a, (mS/m) 


<t„ (mS/m) 


model 19.0 
estimate 24.2 


49.0 
45.5 


200.0 
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Figure 


<r, M (mS/m) 


a, t (mS/m) 


ir„ (in.) 


8 


56.0 (3.0)» 


3.0(1.0) 


0.18(2.0) 


9 


5.0 (10.0) 


2.0 (0.4) 


0.17(1.0) 


10 


10.0 (5.0) 


4.0 (0.2) 


0.21 (0.6) 


11 


11.0 (8.0) 


0.2 (0) 


0.08 (0.4) 


12 


58.0 (2.0) 


140.0(12.0) 


1.2 (4.0) 


13 


3.0 (0.6) 


1.6 (3.0) 


0.18(0.9) 


14 


10.0 (2.0) 


6.7 (14.0) 


1.1 (3.0) 


15 


10.0(0.5) 


1.8 (0.9) 


0.13(0.4) 



Table 1: Standard deviations of the four profile parameters 
From noise study. The additive noise is normally distributed 
noise with 10, 5, 2, and 1 mS/m standard deviation for, re- 
spectively, the 10-, 20-, 30-, and 60-in. curves. 

u M (unitlcss) 

2.3 (8.0) 
0.6 (14.0) 
0.02(1.0) 
1.2 (10.0) 

2.4 (30.0) 
6.1 (60.0) 
0.5 (14.0) 
0.8 (27.0) 

' Numbers in parentheses are the ratio of the sample standard deviation 
estimate to actual value in percent Results are based upon 100 ex- 
periments. 



and the N n resulting sets of parameters are saved. Sample 
statistical estimates of the resulting standard deviations 
for each of the four parameters are computed. The results 
of this statistical noise study for each of the examples in 
Figures 8- 1 5 are summarized in Table 1 . Each entry con- 
sists of the sample standard-deviation estimate, and in 
parentheses to the immediate right is the ratio of the 
sample standard-deviation estimate to actual value in 
percent (computed percent error results rounded to near- 
est 0.1%). Table 1 supports the conclusion that the in- 
version is quite robust for <r, and w. As expected, the 
statistical variations of the shape parameters N and c M 
are larger. N n is 100. ^ a . 

Figures 8-16 illustrate typical accuracy of DCLM in 
known thick-bed environments. Figure 1 7 is a 1 00-ft field 
log recorded with the AIT. Each end of the interval is 
shale, and the center bed is an invaded water-saturated 
sandstone. Note that the DCLM -estimated R t almost 
overlays the 60- and 90-in. curves, which is consistent 
with the rather shallow transition radius w that averages 
around 20 in. The R xo estimate follows the 10-in. log but 
is more resistive and shows more activity. Based upon 
the quality indicators of condition number and the x 2 fit, 
most of these narrow, more resistive R xo zones are prob- 
ably real. In most of this interval, the estimated apparent 
conductivity curves, ^(m), fit the data to several sig- 
nificant figures and have condition numbers, c N , less than 
10' This increase confidence in the results. Two spikes 
at depths of 2,398 and 2,403 ft on the R M curve are 
probably not real because the c N there are in excess of 
10 6 . The ability of the parametric method to predict more 
thin resistive structure in R M than the shallow 10-in. 
curve results from the model-based extrapolation of the 
shape of the deeper transition region towards the well 
bore. 
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Figure 16. Step-invaded center bed of 20- ft thickness. Forward model is based on semianalytic method from Chew et al. (1984). 
In table, bed 1 is uninvaded where z > 10, and bed 2 is center bed where — 10 < z < 10. 



r, (in.) r 2 (in.) <r, (mS/m) a M (mS/m) 



exact, bed 1 


12.0 


12.0 


500.0 


1,000.0 


typical estimate 


9.1 


13.3 


499.3 


1,077.6 


exact, bed 2 


20.0 


20.0 


100.0 


200.0 


center estimate 


8.4 


16.7 


105.0 


222.0 


exact, bed 3 


0.0 


0.0 


1,000.0 
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CONCLUSIONS AND EXTENSIONS 

This paper utilizes a parametric definition of a smooth 
radial invasion profile with a minimum number of pa- 
rameters defining flushed, transition, and uninvaded 
regions. The parametric profiles are maximally flat in the 
flushed zone and in the limit of large N %o over contin- 
uously to the traditional three-parameter <s xoy <r„ and D, 
profile upon which the tornado-chart interpretation is 
based. Numerical results show that for shallow invasion, 
the accuracy of the a t estimate is better than o xo . The 
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maximally flat-parameter model improves the sensitivity 
of the measurement to tr t . The parameters of the profile 
are shown by examples to be accurately estimated from 
median-log data. A new nonlinear algorithm (DCLM), 
coupled with the maximally flat-parametric invasion 
model, is key to extending the ability to deliver reliable 
R t estimates in difficult invasion cases. The theory is de- 
veloped in sufficient generality to estimate six or more 
parameters where annuli could in principle be detected. 
Induction data alone however are limited to resolving no 
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are radii of invasion r, «^r££££}£Z «Stivity). 
then the 30, 20, 10, and 

more than perhaps si, parameters. Fu^ wo^ js p^ed 
to Primate more general nonmonotonic pronles ana u> 

»j ~£ 
provides sood «. estimates in the presence oran annul 
profile. 

NOMENCLATURE 

The following lists mathematical symbols, short ^defr- 

nitils. and SI units. ^££&:Z5Z£ 
results. I f no unit is provided, the Quantity 
matrix whose units are elemen t depe ndent. 

condition number of L matrix 

model constraint for xH-*), unitless 
diameter of invasion, m 
1/[1 + {p/^n - Butterworth profile, unit- 



kHz = 

L > 
L' w ■ 
m : 

ntnc- 

ntok 

m, 

m, 
WI4 
N 



N p - 
/V/, ■ 

n. - 

R, 
r» 

u(p — Po) 



C.v 
C(ni) 

D, 

Jlp, w. AO 



— w 



H 



- step response finite element matrix., unitless 
= Hessian matrix in quadratic fit to x 
= induction-array index, umUess 
= radial pseudo-response function, ohm 
= radial cumulative pseudo-response func- 
tion, ohm-m 



z — 
P* - 

01 - 

5(p - Po) ~ 

6 m = 
A„ = 
, Act = 

y = 



transmitter frequency in k-° Hc «?; 
normalized radial response function, l/m 
LM stabilized H matrix 
LM matrix in terms of q, parameters 
mode) parameter vector with elements m, 
iteration update model vector m 

IS£L"tl*rSi mode, parnm^r 
nominal value 
<r = estimate of o„, S/m 
; 5, = estimate of c„ S/m 
^ - estimate of w. m 

AT- estimate of slope parameter AT. unittess 
four-parameter slope parameter, umtiess 
num£r of radial points used to determine 

^otStion ^ pseudo-response 
function, ohm-m 

number of model parameters, unitless 
number of radial pixels 
finite-element matrix, unitless 
!j/f^J nonlinear change of variable to en- 
force positivity 
= mud-filtrate resistivity, ohm-m 
= connate-water resistivity, ohm-m 
= formation resistivity, ohm-m 
= flushed-zone resistivity, ohm-m 
= w(l - 2/A0 = radius defining end of flushed 

= ™\'-?2/N) = radius defining beginning of 

undisturbed zone, m 
« step function = 1 for p > Po. and zero oth . .. 

crwisc • 
four-parameter transition midpoint, m 
radiarweight-function filter coefficients, 

^parent-conductivity focusing-filter 
weights, unitless 
borehole axial coordinate, m 
Sear-term coefficient in quadratic approx- 

SXSSb of LM matrix equation in 
terms of ^parameters 
impulse response function, l/m 
Krnocker delta - 1 for n « m. and 0 oth- 
erwise, unitless 

mdr/step change in ^nductivity S/m 
constant in quadratic ^approximation, also 
dvnamic constraint ratio, unitless 
LM s^bilizer and switch parameter, unit- 

less 

triangular finite^lement centered at p - p m . 
unitless 
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«0 = 
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*a\ = 
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o a/ (m) = 
X 2 (m) = 
V « 



<p = azimuthal angle in borehole coordinates* 

radians 
r - radial coordinate, m 
>,„ = radial coordinate at wth pixel, m 
<r,- - standard deviation of noise in ith log cr 0/J 
S/m 

= radial step conductivity profile, S/m 
s dynamic-constraint factor determined by 

equation (45) 
= four-parameter radial conductivity profile, 
S/m 

= six-parameter radial conductivity profile, 
S/m 

: flushed-zone conductivity, S/m 
■ annulus conductivity, S/m 
ith median-log apparent conductivity, S/m 
radial profile conductivity at radius p = p m , 
S/m 

constraint function kth parameter tolerance 
a constant or background conductivity, S/m 
jth raw induction array apparent conduc- 
tivity, S/m 

yth step-response array apparent conductiv- 
ity, S/m 

yth step-response array apparent conductiv- 
ity, S/m 

jth induction array apparent conductivity 
for formation <r 0 , S/m 
model estimate of log o ai> S/m 
formation conductivity, S/m 
data fit parameter, unitless 
xKm) - C(m) 
gradient vector operator 
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APPENDIX 

NORMALIZATION OF 
PSEUDO-RESPONSE FUNCTIONS 

For completeness, the homogeneous background gain 
function JV(<7q) for apparent conductivity as used in equa- 
tion (8) is determined. In the most simple case of a mag- 
netic dipole with coaxial coupling (coils are parallel, cen- 
tered on the z axis, and in planes perpendicular to z axis), 
the received voltage v is from elementary theory (for a 
time dependence e"™) 



where 



v =* k(\ - ikr)^. 



k = 2icjpL 0 m R m r I ( /r 3 , 



(A.1) 



r=\z r ~ z*i, 
m T = vp\N T% 



March- April, 1992 



The Log Analyst 



109 



26/27 * RCVD AT 3/22/2003 4:15:14 PM [Eastern Standard Time) - SVR:USPTO-EFXRF-1/5 ■ DNJS:8729308 ■ CS(D:71 32388008 * DURATION (mnvss): 19-08 



qa/22/2005 15:35 FAI 7132388008 CONLEY, ROSE 11027 



Howard 



Ho = 4xl0- 7 H/m, and 
k = (iwpoffo)*- 

The distances z and r are in cylindrical coordinates, N is 
number of turns, « is the circular frequency, T and R 
subscr-pts denote transmitter and receiver, «r 0 is the back- 
ground conductivity in S/m, and is the permeability 
in H/m. To convert from voltage to apparent conductiv- 
ity the voltage v is divided by the sensitivity 3v/*r|_ 0 . 
where if v is computed numerically the evaluation is made 
at some finite but small conductivity where the relation 
between voltage and conductivity is linear, but the com- 
putation is accurate. Thus the dipole normalization is 



(A.2) 



Use equation (A. 1) to evaluate.equation (A.2) resulting in 
NM,) - < A - 3 > 
If the finite size of the loop antennas is important, the 
more general expression for the receiver voltage in the 
same geometry as equation (A.1) is 

and 



N(e 0 ) 



O 



(1-2 sin 2 0y* dO 



Ic £ (l - 2sin 1 0)i<rf0 



(A. 5) 



where 



0 - Vk* - A 2 , 

Im ft a 0, 

L " \Zfi- z r \, 

k = - I^HLom R mjJ{Trp R Pi), 

and J t is the Bessel function of the first kind of order 1. 
In the limit of small p r and Pr, equation (A.4) becomes 
equal to equation (A. 1). The resulting normalization func- 
tion In this case is much more complicated. Nonetheless, 
it can be shown that in this more general case the nor- 
malization is exactly given by 



R = V^TpVl - 7 J s»n 2 «. 
a - O + p\ + p\, 

6 - 2p«Pr. 

7 = 2p7(« + P)- 

In the limit of small loops, equation (A.5) goes to the 
dipole result (A.3). Multi-coil sonde normalizations are 
simple linear combinations of these results. 
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